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Abstract— This paper presents the time-step and grid-size calculation 
methods for solving time series wave differential equations using the 
Finite Difference Method. 


The methods are formulated using the Taylor series, assuming that at a 
very small time-step and grid-size, the number of terms in the Taylor 
series starting with 2nd order with the highest order is much smaller than 
the order 1. 





I. INTRODUCTION 

The formulation of basic hydrodynamic equations, 
both the continuity equation and the momentum equation, 
is developed using the Taylor series to one derivative 
(Dean (1991)). It is argued that in which the time- 
stepdtand grid size (ôx, Z) is very small, calculation on 
the Taylor series can be done up to one order0 (ô+). 
The solution of the governing equations numerically is by 
using the discrete method, both the Finite Difference 
Method (FDM) and the Finite Element Method (FEM) by 
dividing the domains into grid-points. The grid-point 
formation requires a small grid-size under the formulation 
conditions, assuming the time-stepdtand grid size 
(ôx, 6z)are very small. 
In FDM, the equations of FDM are formulated by 
truncating the Taylor series, using only the first or second- 
order (Smith, G.D (1985)). Therefore, the application of 
the equations isrequired to use a very small time-step and 
grid-size where the Taylor series can be used in order 1 
and order 2 only. 
Thus, it takes a small time-step and grid-size following the 
formulation of governing equations as well as in the 
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formulation of the FDM equation in which FDM is used in 
the solving governing-equations. 

The accuracy of the Taylor series is determined not only 
by the length of the terms used but also by the grid-size 
used, where the smaller the grid size, the better the 
accuracy. At very small time-steps and grid-sizes, the 
number of terms of order 2 and greater is much smaller 
than terms of order 1. This, this condition is the basis for 
the development of the grid-size calculation method in this 
study. 

The grid-size study is carried out using the sinusoidal 
function, where the water wave equationsare the result of 
an analytical solution to the Laplace’s equation by using 
the variable separation method is a sinusoidal function in 
time and space. 


Il. THE SETTING OF TIME STEPStAND GRID 
SIZE 6x 

The water wave equations are obtained analytically, 
by solving Laplace's equation with the variable separation 
method in the form of multiplication between the 
sinusoidal function and the hyperbolic function (Dean 
(1991)). The horizontal axis-x and time-t form a 
sinusoidal function while the vertical axis-z form a 
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hyperbolic function. By substituting the water wave 
equations into the Taylor series, the equations for the time- 
step p Stand grid-size dxand6éz were obtained. 

The formulation is carried out by two methods. The first 
method isa separate method where the time-stepdtis 
calculated first and then inputdét, grid-size dxandézare 
calculated. The second method is to substitute the relation 
betweendx, 6z, and the time step dtin the Taylor series to 
obtain the equation fordt. Then, dxandézarecalculated 
using the relation betweendxanddézwith the obtained ôt. 


2.1. Separate Calculation 

a. ôt Calculation 

The Taylor series form (Thomas, George B., Finney, Ross 
L.(1996)) for a function of time-t is, 


ôt) = ô a 
f(tt+6t) =f) + tae 


ôt? d*f  6t2 d?f | S&t* a*f 
2 dt? 6 dt? 24 dt* 
. 1) 
To enable Taylor series can be used up to one derivative 
only, then 
ôt? ad?f | St2d?f | St*da*f 
2 dt? 6 dt? 24 dtt 
d 
st” 
dt 





<e 





ôtin numerator and denominator cancel each other out, 

dea? f St d?f StFatf 

2 dt? 6 dt? 24 dtt 
df 














dt 
Used 
f(t) = Acos ot.... (3) 
o= TT, Angular frequency, T: period 
Substitution (3) to (2) where the amplitude A on the 
numerator and denominator cancel each other out is done 
at the characteristic point and is done at the characteristic 
point, cos ot = sin ot, in which the terms 
cos otandsin otin the numerator and denominator also 
cancel each other out. 

ôt 00" 4 Ot 4 
- rad = gl + 34° Se 
Since the term in the absolute sign | | is negative, then 
the positive value of the equation becomes, 
êt ôt? ôt? 





2° + 6 o? = g” <e 
By taking the equal sign, 
2 3 
ee — ËE o = gn (4) 
2 6 24 


Equation (4) can be solved by the Newton-Raphson 
method (Allen, Myron B.; Isaacson, Eli L.(1998)). The 
results of the calculation of tusing equation (4), for 
several wave periods, are presented in Table (2.1), with an 
accuracy level ofe = 0.001. Meanwhile, Table (2.2) 
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presents the calculation result of ôt for the wave period is 
8 sec. with several levels of accuracy. 
Table.2.1: The value of ôt for severalT values withe 0.001 












































T ôt bt 
(sec.) (sec.) T 
6 0,00191 | 0,00032 
7 0,00223 | 0,00032 
8 0,00254 | 0,00032 
9 0,00286 | 0,00032 
10 0,00318 | 0,00032 
11 0,0035 | 0,00032 
12 0,00382 | 0,00032 
13 0,00414 | 0,00032 
14 0,00445 | 0,00032 
15 0,00477 | 0,00032 








Table.2.2: The value of Stat several values ofe, T = 8 sec. 

















E€ ôt bt 
(sec.) T 
0,1 0,23899 | 0,02987 
0,01 0,0253 | 0,00316 
0,001 | 0,00254 | 0,00032 
0,0001 | 0,00025 | 0,00003 

















Table (2.1) shows that Z is constant for all periods of T for 


a given level of accuracye. Meanwhile, Table (2.2) shows 
the change in ôt along with the change in the level of 
accuracy€. 


b. 6x Calculation 

The calculation of ôx on the function f(x, t)is done using 
ôt as the input obtained by the calculation method as 
discussed in the previous section. In this section, the 
Taylor Series is only used up to the 3rd differential 
or0 (6%), considering that thedtvalue obtained is quite 
accurate, and to facilitate writing. Taylor’s series up to the 
third derivative for a function f (x, t)is, 


f(x + 6x,t + 6t) = f(x,t) + s1 + 52 + 53..... (5) 
Where, 

= 5a ar 
Sı = ôt aE + 6x prone (6) 

7 ay, aer 
~ 2 at? atdx 2 ax? 
ôt? adf = St? d3f 
= —— + 6x —— 

6 dt? 2 dt*dx 
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ôx? adf dx d3f 
+6t ———~ + — 
2 dtdx? 6 dx? 


.. (8) 


To enable the Taylor series used with only one derivative, 


then 
S2+S3 








S41 


The function below is used, 
f (x,t) = Acos kx cosot..... (10) 


By using the relation (6) and (7), 
ot? 8? f StS af — ôx? af 





Sy st + x 
dt x 


Substitutef (x,t) done at a characteristic point 





wheresinot = cosot = sinkx = coskx, 
2 2 
So oe 0 Stéxok +e k? (11) 


Sy ôto+ôxk 
In the same way, it was obtained 


6t3 bt? 8x? 8x? 
S3 o? 5x02 k-6t?*ok?-2* k3 
—_6 2 2 6 (12) 


Sy ôto+ôxk 








Intuitively, it can be estimated that (11) is positive, while 
(12) is negative. However, considering that (11) is greater 
than (12), the sum of the two equations will be positive. 


Substituting the two equations to (9) and assuming that the 
term in the absolute sign is positive, the equation is 
obtained. 


Se ôt? 2, ot 5 
Cee gee ere 


ôt? 
= Sto +- a’ +e kôx 


k? k? 
+(1 + dto) = dx" = |G oe = 


(13) 


ôxcan be calculated by (13) with ôt as the input, 
wheredtis obtained from the previous procedure. 


Table (2.3) presents the calculation results of dx for 
several T wave periods, at a water depth of 10m. The 
wavelength is calculated using the dispersion equation 
from the linear wave theory according to Dean (1991).The 
level of accuracy used ise = 0.001. 
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Table.2.3: 5x Calculation 
T ôt ôx L 





ôx 

(sec.) (sec.) (m) (m) L 
6 0,00191 | 0,04643 | 48,4062 | 0,00096 
7 0,00223 | 0,05738 | 59,8212 | 0,00096 
8 0,00254 | 0,068 | 70,8984 | 0,00096 
9 0,00286 | 0,07839 | 81,7267 | 0,00096 
10 | 0,00318 | 0,0886 | 92,3739 | 0,00096 
11 0,0035 | 0,09868 | 102,887 | 0,00096 
12 | 0,00382 | 0,10867 | 113,299 | 0,00096 
13 0,00414 | 0,11858 | 123,633 | 0,00096 
14 | 0,00445 | 0,12843 | 133,905 | 0,00096 
15 0,00477 | 0,13824 | 144,128 | 0,00096 





















































The calculation results show that the value of dxobtained 
is quite small compared to the wavelength L, which shows 
that the value is quite relevant to use. Since the value 


o — might be too small, the comparison between the wave 
phase velocity is reviewed, C = =, and-=(Table (2.4)). 


Table.2.4: Comparison betweenandC 



































ôx 
T ôt C orst 
(sec.) (m/sec.) | (m/sec) C 
6 24,3257 | 8,0677 | 3,01519 
7 25,7675 | 8,54589 | 3,01519 
8 26,7215 | 8,86229 | 3,01519 
9 27,3802 | 9,08074 | 3,01519 
10 27,8525 | 9,23739 | 3,01519 
11 28,2022 | 9,35337 | 3,01519 
12 28,4682 | 9,44158 | 3,01519 
13 28,6751 | 9,51022 | 3,01519 
14 28,8392 | 9,56465 | 3,01519 
15 28,9716 | 9,60854 | 3,01519 




















Table (2.4) shows that—is much greater thanCthere the 


ôx 
ratio t = 3.015190rdx = 3.01519 C ôt. This is 


following with the Courant criteria,ĝx = 3.0 C ôt 
(Courant (2928)). 
2.2. dtcalculation by substitutingdx. 
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In this section dt and dxare processed using the same 
procedure as the previous one where 6dxis substituted 
withdt, in which 

ôx = yCét = yzôt Ers (14) 


C= nis wave celerity 
L= Wavelength 


y:Coefficient (when used the results in Table (2.4.),y = 
3.01519) 


Furthermore, using the relation (6) and (7), wheredxis 
substituted by (14) to obtain: 

2 
: = (=) oot ......(15) 
By using the relations (6) and (8) and by substitutingédx 
with (14), 


3 


2 
Sa (“eer o26t?.... (16) 





Sy 1+y 


Substitute (15) and (16) into (9) and taking a value equal 
to, 


1 y,v,¥ 
eas a 


E | 076? 
1+y ° 
1y Ê 
1+y 


By inputtingo, candy, dtcan be calculated with (14). 
Meanwhile, y can be obtained using the previous 
calculation results,y = 3.01519. 


The calculation of ôx in the table (2.5) is done in the water 
depth h = 10 m, while the wave number kis calculated 
using the dispersion equation of the linear wave theory 
(Dean (1991)). It shows that ôt and dx in Table (2.5) are 
more or less the same as the calculation results in Table 
(2.3). However, the calculation of ôt and (17) was 
preferable considering that there is a clear interaction 
between dt andéx. 


Table.2.5: The results of the calculation of ôt and 6x with 
(17), fory = 3.01519. 

















T ôt ôx 5 
x 
(sec) | (sec) (m) T 
6 0,00191 | 0,04619 | 0,00095 
7 0,00223 | 0,05709 | 0,00095 
8 0,00255 | 0,06766 | 0,00095 
9 0,00286 | 0,07799 | 0,00095 
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10 0,00318 | 0,08815 | 0,00095 
11 0,0035 | 0,09819 | 0,00095 
12 0,00382 | 0,10812 | 0,00095 
13 0,00414 | 0,11798 | 0,00095 
14 0,00445 | 0,12779 | 0,00095 
15 0,00477 | 0,13754 | 0,00095 








2.3. 6z calculation in the functionf (x, z, t) 


This section formulates the calculation method 6z for a 
function f (x, z,t). The calculation method is as discussed 
in the previous section, where dt and dxare the inputs 
obtained from the previous procedure on a sinusoidal 
function f (x,t). 


Taylor series for a function f (x, z, t) is, 
f(x + ôx,z + ôz,t + ôt) = 
f(%,z,t) + s1 +52 +53 
we (19) 


sı = ôt £4 6x44 52%... 20) 
dt dx dz 











ôt? a2 f  éx2a2 
ee Sa gi ay 
2 de dxdt 2 dx? 
af af szaf 
ôt 5x6. pelea 
totoz ge a a az? 
ee (21) 


sear St? af 
S3 = + x 
6 dt? 2 0” dxae2 
st? g 8x2? a 
Oey aac ee a 
2° azde? 2 dxdt 
ôz? af af 
t ôtôxő 
2 azat  0O~OZ azat 
8x3 If 8z? a? 
gee Af pora 
6 dx? 6 az? 











ôzwill be calculated using equation (9). 


As a functionf (x, z, t), the following form of function is 
used 


f(x, z,t) = coshk(h + z)coskxcosot.... (23) 


At the characteristic point, the sinusoidal terms in the 
numerator and denominator might cancel each other out. 
The equation is done atz = 7, whereyis water surface 
elevation and is done in deep water where tanhk(h + 
n) = lorsinhk(h +n) = coshk(h +7). Thus, the 
hyperbolic function of the numerator and the denominators 
might eliminate each other.By using the wave-number 


Page | 283 


Madson Rodrigues Lemos et al. 


conservation law (Hutahaean (2020)), the 


relationtank(h + 7) = 1remains valid in shallow waters. 


After removing the sinusoidal and hyperbolic elements, the 
equations of S4, S2 and sz are, 


Sı = —6to — xk + 62k... (24) 


2 6x2 


Ob 24 5 
so = RES + ôtôxok gt 


6 2 
—Stőzok — Sxőzk? + Lg 


ae (25) 
ot a JOU 2% stt 2% 
53 = 7 z Oxo z ZO 
pe ge 
27 a? 
bz? 
—dt— ok? + 6tdxdzok? 
ôx’ 623 
— k? + — k?’ 
me Re 
oe (26) 





Substitutings,, S2, ands,to (9) assuming thats 
1 
positive and using the equal sign, the equation fordzis in 

the form of a 3-degree polynomial. 

The results of the calculation of 6z where ôt and xare 
calculated using the procedure in sub-chapter (2.2) are 
presented in Table (2.6), where the calculation at a water 


depth of 10m using the accuracy levele = 0.001. 
Table.2.6: 6z Calculation 





T ôt ôx Oz 



































Oz 

(sec) (sec) (m) (m) ôx 
6 0,00191 | 0,04619 | 0,13859 | 3,00012 
7 0,00223 | 0,05709 | 0,17127 | 3,00012 
8 0,00255 | 0,06766 | 0,20298 | 3,00012 
9 0,00286 | 0,07799 | 0,23398 | 3,00012 
10 0,00318 | 0,08815 | 0,26447 | 3,00012 
11 0,0035 | 0,09819 | 0,29457 | 3,00012 
12 0,00382 | 0,10812 | 0,32438 | 3,00012 
13 0,00414 | 0,11798 | 0,35396 | 3,00012 
14 0,00445 | 0,12779 | 0,38337 | 3,00012 
15 0,00477 | 0,13754 | 0,41264 | 3,00012 
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Something is interesting enough to note that the value of 
6 : : P 
i = 3.0001 2isconstat for allwave periods. Thus, it can be 


assumed that: 
2 
6z = wa = y?C6t.... (27) 


Where y = 3 can be used. 


Hence, a calculation method might be formulated where 
ôx and dzaresubstituted by dt using the relations in (14) 
and (27) to (24), (25), and (26). 


sı = (—1 — y + y? )oôt... (28) 


1 3 1 
sz = (4+ y = žy? -7° +274) 0761... (29) 





S3 7y3 -ży* +1y6) a?ôt?.... (30) 


Ga 


Substituting (28), (29), and (30) to (9) and assuming that 
the term in the absolute sign is positive, the quadratic 
equation of dtis obtained as follows, 


a 8t? + b ôt + c = 0...... (31) 





a = (+ ży 423-274 + 2°) 2... (32) 
pi 3 1 

b=(-S+y-27?-73 +14) 0... (33) 

c=-(-1-yty?)e... 34) 


The results of the calculation of ôt with (31) for several 
wave periods, at a water depth of 10m with an accuracy 
level ofe = 0.001, whiley = 3.00012are presented in 
Table (2.7). 


Table.2.7: Calculation results of ôt, 6x and 6z with (31) 















































T ôt ôx Oz 
(sec) (sec) (m) (m) 
6 0,00176 | 0,04261 | 0,12782 
7 0,00205 | 0,05265 | 0,15797 
8 0,00235 | 0,0624 | 0,18722 
9 0,00264 | 0,07193 | 0,21581 
10 0,00293 | 0,08131 | 0,24392 
11 0,00323 | 0,09056 | 0,27169 
12 0,00352 | 0,09972 | 0,29918 
13 0,00381 | 0,10882 | 0,32647 
14 0,00411 | 0,11786 | 0,35359 
15 0,0044 | 0,12686 | 0,38059 








There is a relatively small difference between the results in 
Table (2.7) and the results in Table (2.6). It can be said 
that with (31) simultaneous calculations onôt, dxandéz, 
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the results of calculation with (31), in Table (2.7) will be 
more reliable than the results of calculations in (2.6). 


ee 6 6 : 
It is important to note that and correlate with the 


velocity of the C wave phase, in equations (14) and (27), 
in interpreting these two parameters. 


MI. AIRY’S LONG WAVE EQUATION 
COMPLETION 

In this section, Airy’s Long Wave Equation or also 
known as the Shallow Water Equation is resolved. The 
equationis actually for small amplitude and longwave. To 
be used for short and high amplitude waves, a modification 
of the equation might be used. The formulation is not 
presented in this study considering that this study does not 
aim to develop a water wave model, but only to test the 
reliability of the time-step and grid-size produced in the 
previous sections. 
Governing equations consists of two equations, namely 


a. Water surface equation 

a 1 aUh a 

2 ee E025 2... (35) 
dt 0.5m dx dx 

b. Momentum equation 

au 1 /3auu a 

Z E + g2)... 66) 

at 9 \2 ax ax 


7: water surface elevation 

U: depth average particle velocity 
h: still water depth 

g: gravitation acceleration 





n(x, t) 








Oooo 


Fig.1: The axes system and variables in SWE 


Both equations are solved by modifying MacCormack’s 
predictor-corrector method (Anderson, J.D. Jr.(1994)). 
The original predictor-connector method is a combination 
of the forward difference with the backward difference in 
solving the time differential (Anderson, J.D. Jr.(1994)). In 
this study, a combination of central and backward 
differences was used. 

The simulation results for waves with a period of 9 sec., at 
constant still water depth h =10 m, wave amplitude 0.8 m 
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on execution for 10 times the wave period or 80 sec., are 
presented in Fig. (2). Time-stepôtand grid-size xare 
generated using the accuracy level ofe = 0.001. 


1 
0.5 
£ o 

low 

-0.5 
-1 

0 100 200 300 

x (m) 
eta eta-max 








Fig.2: Model results in constant water depth,h = 10 m. 


Fig. (2) shows that the wave curve is stable at execution 
for 80 sec. Thus, the wave crest elevation Nmax Was stable. 
This shows the stability of the model, the numerical 
method used, as well as the good time-step and grid-size. 
In case the time-step and grid-size are not correct, for 
example, the grid-size is too large, and then a reduction in 
waveamplitude or wave height might occur. 

Given the wavelength of the model is shorter than the 
wavelength of the linear wave theory, the calculation of 
the grid-size ôx used the following dispersion equation to 
obtain the appropriate grid-size: 

o? = gk tanh(kh)....... (37) 

However, this study does not propose a new dispersion 
equation, but only the calculation of the grid-size. 
Furthermore, the model is worked on the sloping bottom 
where the initial still water depth is 10m, while the final 


still water depth is 1m. The bottom slope is = with a 


channel length of 300m. The wave period is 8 sec, with a 

wave amplitude of 0.8 m. The model results are presented 
in Fig (3). 
1:5 
1 
_ 05 
E o 
-0.5 
-1 
-1.5 

0 100 200 300 


x (m) 








eta eta-max 


Fig.3: Results of model execution at the sloping bottom. 
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Fig (3) visualizes the occurrence of shoaling. Meanwhile, 
in the very shallow waters, there is a large increase in 
wave crest which possibly breaking the model. The 
concern, in this case, is the stability of the model in 
unstable wave conditions, at the time of breaking, which 
also shows that the time-step and grid-size used are quite 
good. 


IV. CONCLUSION 

By using the right time-step and grid-size, the 
Taylor series can be cut using only one derivative (order 
1). Thus, the equations in FDM are exact. 
Moreover, the use of the right time-size and grid-size, the 
solution of the governing equations of the hydrodynamic 
equations is under the conditions of the formulation, where 
generally the governing-equations are formulated by 
cutting the Taylor series in order 1 only. 
It is important to note that the results of this study are that 
the grid-size onthe Taylor series in a sinusoidal function is 
correlated with the phase velocity instead of the particle 
velocity or the current velocity. This should be considered 
in formulating the governing equations of the wave model, 
especially in the formulation of the momentum equation. 
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